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Abstract 

We present a simple construction method for Feller processes and a framework for the generation of 
sample paths of Feller processes. The construction is based on state space dependent mixing of Levy 
processes. 

Brownian Motion is one of the most frequently used continuous time Markov processes in applications. 
In recent years also Levy processes, of which Brownian Motion is a special case, have become increasingly 
popular. 

Levy processes are spatially homogeneous, but empirical data often suggest the use of spatially inho- 
mogencous processes. Thus it seems necessary to go to the next level of generalization: Feller processes. 
These include Levy processes and in particular Brownian motion as special cases but allow spatial inho- 
mogencities. 

Many properties of Feller processes are known, but proving the very existence is, in general, very 
technical. Moreover, an applicable framework for the generation of sample paths of a Feller process was 
missing. We explain, with practitioners in mind, how to overcome both of these obstacles. In particular 
our simulation technique allows to apply Monte Carlo methods to Feller processes. 

Introduction 

The paper is written especially for practitioners and applied scientists. It is based on two recent papers 
in stochastic analysis [TJOH]. We will start with a survey of applications of Feller processes. Thereafter 
we recall some existence and approximation results. In the last part of the introduction we give the 
necessary definitions. 

The main part of the paper contains a simple existence result for Feller processes and a description 
of the general simulation scheme. These results will be followed by several examples. 

The source code for the simulations can be found as supporting information (Appendix SI). 

Motivation 

Brownian motion and more general Levy processes are used as models in many areas: For example in 
medicine to model the spreading of diseases [23], in genetics in connection with the maximal segmental 
score [13] . in biology for the movement patterns of various animals (cf. |31| and the references therein), 
for various phenomena in physics |37j and in financial mathematics |10j . In these models the spatial 
homogeneity is often assumed for simplicity, but empirical data or theoretical considerations suggest that 
the underlying process is actually state space dependent. Thus Feller processes would serve as more 
realistic models. We give some explicit examples: 
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• In hydrology stable processes are used as models for the movement of particles in contaminated 
ground water. It has been shown that state space dependent models provide a better fit to empirical 
data [3nil3H] • Also on an intuitive level it seems natural that different kinds of soils have different 
properties. Thus the movement of a particle should depend on its current position, i.e. the soil it 
is currently in. 

• In geology also stable processes are used in models for the temperature change. Based on ice-core 
data the temperatures in the last-glacial and Holocene periods are recorded. Statistical analysis 
showed that the temperature change in the last-glacial periods is stable with index 1.75 and in the 
Holocene periods it is Gaussian, i.e. stable with index 2 (see Fig. 4 in [12]). 

• For a technical example from physics note that the fluctuations of the ion saturation current mea- 
sured by Langmuir probes in the edge plasma of the Uragan-3M stellarator-torsatron are alpha- 
stable and the alpha depends on the distance from the plasma boundary |15) . 

• Anomalous diffusive behavior has been observed in various physical systems and a standard model 
for this behavior are continuous time random walks (CTRWs) [27j|33]. To study these systems the 
limiting particle distribution is a major tool, which is in fact a Feller process pQ. 

• In mathematical finance the idea of extending Levy processes to Levy-like Feller processes was 
first introduced in [2]. The proposed procedure is simple: A given Levy model usually uses a 
parameter dependent class of Levy processes. Now one makes the parameters of the Levy process 
(in its characteristic exponent) price-dependent, i.e. the increment of the process shall depend on 
the current price. This procedure is applicable to every class of Levy processes, but the existence 
has to be shown for each class separately [2]|3l|6]. 

Thus there is plenty of evidence that Feller processes can be used as suitable models for real-world 
phenomena. 

Existence and Approximation 

Up to now general Feller processes were not very popular in applications. This might be due to the fact 
that the existence and construction of Feller processes is a major problem. There are many approaches: 
Using the Hille-Yosida theorem and Kolmogorov's construction p~8j[20], solving the associated evolution 
equation (Kolmogorov's backwards equation) [4l[5][24j[25| , proving the well-posedness of the martingale 
problem [511180 35]. solving a stochastic differential equation [TW2"2"l[3l)] . The conditions for these construc- 
tions are usually quite technical. Nevertheless, let us stress that the proof of the very existence is crucial 
for the use of Feller processes. Some explicit examples to illustrate this will be given at the end of the 
next section. 

Our construction will not yield processes as general as the previous ones, but it will still provide a 
rich class of examples. In fact the presented method is just a simple consequence of a recent result on 
the solutions to certain stochastic differential equations [34] . 

Furthermore each of the above mentioned methods also provides an approximation to the constructed 
Feller process. Most of them are not usable for simulations or work only under technical conditions. Also 
further general approximation schemes exist, for example the Markov chain approximation in [26] . But 
also the latter is not useful for simulations, since the explicit distribution of the increments of the chain 
is unknown. 

In contrast to these we derived in [7] a very general approximation scheme for Feller processes which 
is also usable for simulations. We will present here this method for practitioners. 
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Levy processes and Feller processes 

Within different fields the terms Levy process and Feller process are sometimes used for different objects. 
Thus we will clarify our notion by giving precise definitions and mentioning some of the common uses of 
these terms. 

A stochastic process is a family of random variables indexed by a time parameter t £ [0, oo) on a 
probability space (f2, J- ', P). For simplicity we concentrate on one-dimensional processes. The expectation 
with respect to the measure P will be denoted by E. 

Although this will not appear explicitly in the sequel, a process will always be equipped with its so-called 
natural filtration, which is a formal way of taking into account all the information related to the history 
of the process. Technically the filtration, which is an increasing family of sigma fields indexed by time, 
is important since a change from the natural filtration to another filtration might alter the properties of 
the process dramatically. 

A Levy process starting in Lq := is a stochastic process (£t)te[o,oo) with 

- independent increments: The random variables L tl , L t2 — L tl , L t3 — L t2 , . . . are independent for 
every increasing sequence {t n ) nl =m, 

- stationary increments: L t — L s has the same distribution as L t - S for all s < t, 

- cadlag paths: Almost every sample path is a right continuous function with left limits. 

For equivalent definitions and a comprehensive mathematical treatment of Levy processes and their 
properties see [52] . 

Note that the term Levy flight often refers to a process which is a continuous time random walk 
(CTRW) with spatial increments from a one-sided or two-sided stable distribution (the former is also 
called Levy distribution) . In our notion the processes associated with these increments are Levy processes 
which are called stable subordinator and stable process, respectively. 

A Levy process (£t)te[o,oo) on its probability space is completely characterized by its Levy exponent 
£ ^ calculated via the characteristic function 



almost every sample path is continuous. In general, Levy processes have discontinuous sample paths, some 
examples with their corresponding exponents are the Poisson process (^(£) = e 1 ' — 1), the symmetric te- 
stable process = with a £ (0, 2]), the Gamma process (ip(£) = ln(l-H£)) and the normal inverse 



Gaussian process (?/>(£) = -ifj£ + 5[^a 2 - {fi + iQ 2 - y 'a 2 - P 2 ) with < \0\ < a, 5 > 0, /i£ R). 



Classes of Levy exponents depend, especially in modeling, on some parameters. Thus one can easily 
construct a family of Levy processes by replacing these parameters by state space dependent functions. 
Another approach to construct families of Levy processes is to introduce a state space dependent mixing 
of some given Levy processes. We will elaborate this in the next section. 



{' t Px)x£R> we can construct for fixed xq £ K, T £ [0, oo), n £ N a Markov chain as follows: 

1. The chain starts at time in xq. 

2. The first step is at time i and it is distributed as L\°\ The chain reaches some point x\. 

3. The second step is at time — and it is distributed as L^\ The chain reaches some point x%- 





i.e. given a family of characteristic exponents 
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4. The third step is at time — and it is distributed as L^ 2 \ The chain reaches some point 23. 

5. ... etc. until time T. 

This Markov chain is spatially inhomogcncous since the distribution of the next step always depends 
on the current position. If the chain converges (in distribution for n — > 00 and every fixed T € [0, 00)) 
then the limit is - under very mild conditions (see [5J and also Theorem 2.5 by [H]) - a Feller process. 
Formally, a Feller process is a stochastic process (^t)tg[o,oo) such that the operators 

T t f(x) :=E(f(X t )\X Q =x), t€ [0,oo),x &R 

satisfy 

Tq = id, T s oT t — T s+t (s,t > 0) (semigroup property) 



and 



lim sup \T t f{x) — f(x)\ = (strong continuity) 



for all / which are continuous and vanish at infinity. 

A Feller processes is sometimes also called: Levy-type process, jump-diffusion, process generated by 
a pseudo-differential operator, process with a Levy generator or process with a Levy-type operator as 
generator. Note that in mathematical finance often the Cox-Ingersoll-Ross process [TTJ is called the Feller 
process, but in our notion this is a Feller diffusion in the sense of |14j . For a comprehensive mathematical 
treatment of Feller processes and their properties see [30] . 

The generator A of a Feller process is defined via 



lim sup 

t— >0 tC i 



Af(x)-^- f M 



t 



= 



for all / such that the limit exists. Moreover, if the limit exists for arbitrarily often differentiable functions 
with compact support then the operator A has on these functions the representation 



Af(x) = - / e** MO / e-^f(y)dy d£, 

Jr Jr 

where for each fixed idl the function £ i— > ip x (£,) is a Levy exponent. Thus a family of Levy processes 
with Levy exponents (ipx)xeR corresponds to the Feller process (^t)te[o,oo) with generator A as above. 

If the corresponding family of Levy processes is a subset of a named class of Levy processes, one calls 
the Feller process also by the name of the class and adds -like or -type to it. Thus for example a Feller 
process corresponding to a class of symmetric stable processes is called symmetric stable-like process. 

In general, as mentioned in the previous section, the construction of a Feller process corresponding to a 
given family of Levy processes is very complicated. It even might be impossible as the following examples 
show: Let (Zj ^)t>o be the family of Levy processes with characteristic exponents ip x (0 = — iffl(ai)£, i.e. 
the Levy processes have deterministic paths = a(x)t. Now if a(x) = x a corresponding Feller process 
exists, starting in x it has the path X t = xe l . But for a{x) = x 2 and a(x) = sgD.(x)^/x a corresponding 
Feller process does not exist, a(x) = x 2 yields paths which do not tend to negative infinity as x — » —00 
and a(x) = sgn(a;) v / x yields paths which are not continuous with respect to the starting position. 

However, we will present in the next section a very simple method to construct Feller processes. 
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Results and Discussion 

Construction of Feller processes by mixing Levy processes 

Suppose we know (for example based on an empirical study) that the process we want to model behaves 

like a Levy process (l| 1 ') 4 >o in a region K\ and like a different Levy process (Lj )t>o in a region K 2 . 
Then we know that a Feller process which models this behavior exists by the following result: 
Theorem. If the sets Ki, K 2 are uniformly separated, i.e. there exists an e > such that 

inf II2; — y\\ > e 

then there exists a Feller process (X t )t>o which behaves like on K\ and like L^ 2 ' on K 2 . 

Proof. Let ip^ be the characteristic exponent of (Li t )t>o for i = 1,2. Under the above condition there 

exist non-negative bounded and Lipschitz continuous functions a^- 1 ' and aS 2 > such that 

a {1 \x) = 1 for all x G K x and a (1) (x) = for all x G if 2, 
a (2) (x) = for all a; G #1 and a {2) {x) = 1 for all x £ K 2 . 
Now set for (i,5,J;£l 

V>( I ) ^^(fO + ^te), *(*) := ( j^J ) T eK lx2 

and note that for a;, £ £ M 

^($ t (.t)0 = 4^ (a {1 \x)£}j + ^ 2 > (aW(x)^ 
holds. Thus corresponding to the family of Levy processes defined by the Levy exponents 
4< x (0 := (x)t) + 4 {2) (a^{x)£) = 

there exists a Feller process as a consequence of Corollary 5.2 from [34] and 4> x (£) = for x £ K, 

holds (i = 1, 2), i.e. (^t)t>o behaves like on Ki for i = 1,2. ■ 

Note that the theorem extends to any finite number of Levy processes (L[ l )t>o (i = l,~,ri) with 
corresponding regions Ki. More generally for any finite number of independent Levy processes (L[ )t>o 
(i = l,..,n) with corresponding characteristic exponents ip^ and non- negative bounded and Lipschitz 
continuous functions x i-> a^'(x) the family (ip x ) X £WL with 

:= (a' 1 ^) + ^ 2 > (a (2) (a:)e) + • • • + ^ (n) (^(a^) 
defines a family of Levy processes ((ifhoo) an d there exists a corresponding Feller process (X t ) t >o- 

V ~ /x€l 

To avoid pathological cases one should assume 0S 1 ' (x) + (x) + . . . + a'"' (x) > for all x. Further 
note that the following equality in distribution holds for all x and t 

if ± (x)L { t 1] + ( x )l[ 2) + ...+ a W ( x )L { t n) . 

Thus if one knows how to simulate increments of the L\ one can also simulate increments of L\ . We 
will see in the next section that simulation of increments of the corresponding family of Levy processes 
is the key to the simulation of the Feller process. 
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Simulation of Feller processes 

Given a Feller process (X t ) te r oo) with corresponding family of Levy processes ( )tefO ool ) we can 
use the following scheme to approximate the sample path of X t : 

1. Select a starting point xo, the time interval [0, T] and the time-step size h. 

2. The first point of the sample path is x = xq at time t = 0. 

3. Draw a random number z from the distribution of (a; is the current position of the sample 
path). 

4. The next point of the sample path is x ~> .t + z at time t t + h. 

5. Repeat 3. and 4. until i 6 (T - /i,T]. 

The simulated path is an approximation of the sample path of the Feller process, in the sense that 
for h — >■ it converges toward the sample path of the Feller process on [0, T]. To be precise, for the 
convergence the Feller process has to be unique for its generator restricted to the test functions and 
the family of Levy processes L\ has to satisfy some mild condition on the x-dependence: The Levy 
exponent ip x (£) has to be bounded by some constant times 1 + |£| 2 uniformly in x, see [7] for further 
details. This condition is satisfied for many common examples of Feller processes, in particular for the 
processes constructed in the previous section. 

The reader familiar with the Euler scheme for Brownian or Levy-driven stochastic differential equa- 
tions (SDEs) will note that the approximation looks like an Euler scheme for an SDE. In fact it is an 
Euler scheme, but the corresponding SDE does not have such a nice form as for example the Levy-driven 
SDEs discussed in [28]. This is due to the fact that in their case for a particular increment all jumps of 
the driving term are transformed in the same manner, but in the general Feller case the transformation 
of each jump can depend explicitly on the jump size. More details on the relation of this scheme to an 
Euler scheme can be found in a forthcoming paper [8]. 

Examples 

We will now present some examples of Feller processes together with simulations of their sample paths. 
The first example will show the generality of the mixture approach, the remaining examples are special 
cases for which the existence has been shown by different techniques. 

All simulations are done with the software package R and the source code of the figures can be 
found as supporting information (Appendix SI). 

Brownian-Poisson-Cauchy-mixture Feller process 

To show the range of possibilities which are covered by the mixture approach we construct a process 
which behaves like 

Brownian motion on (—oo, —6), 

a Poisson process on (—4,4), 

a Cauchy process on (6,co). 

For this we just define a family of Levy processes by the family of characteristic exponents (ip x )x€R with 

^(0 := ai(z)i|£| a + a 2 (x)(l - e"«) + a 3 (a;)|£| 
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where 



ai(x) 



and 



■r+6 
2 



if x < —6 
if iG (-6, 
otherwise 



-4) , a 2 (x) 








if x > 6 
if x e (4,6) 
otherwise 



if xe [-4,4] 
if xG (-6,-4) 
if x e (4, 6) 
otherwise 



These functions are Lipschitz continuous and thus a corresponding Feller process exists. Figure [T] 
shows some samples of this process on [0, 20] with time-step size y^. One can observe that the process 
behaves like a Poisson process around the origin, like a Cauchy process above 6 and like Brownian motion 
below -6. 



Symmetric stable-like process 

A Levy process L t is a symmetric-a-stable process if there exists an a € (0,2] such that its characteristic 
function is given by 

If we now define a function x i— > a{x) where a(x) takes only values in (0, 2] then there exists a family of 
of Levy processes (£^) te [o.oo) such that for fixed x the the process lJf' has the characteristic function 

A corresponding Feller process exists and is unique if the function x i— > a(x) is Lipschitz continuous and 
bounded away from and 2 [3]. 

Figure [2] shows some samples of a stable-like Feller process on [0, 20] with time-step size jq and 



, . -,19 f fx 
a{x) :=1+ lo((4- 



i.e. x i— » a(x) is a function which is Lipschitz continuous (but not smooth) oscillating between 1 and 
(nearly) 2. To understand the figure note that we color coded the state space: red indicates awl, 
yellow indicates a 2 and the values between these extremes are colored with the corresponding shade 
of orange. Now one can observe that the process behaves in the red areas like a Cauchy process and the 
more yellow the state becomes, the more the process behaves like Brownian motion. 



Normal inverse Gaussian-like process 

The characteristic function of a normal inverse Gaussian process L t is given by 
Ee liLt = cxp (tifi - tS (y<* 2 - (P + W - Va 2 - /3 2 )) 

where a > 0, —a < /3 < a, < S and p.. If we replace a, /3, 5, /i by arbitrarily often diffcrcntiable 
bounded functions a(x), P(x) 1 8{x), fi(x) with bounded derivatives and assume that ther exist constants 
c, C > such that S(x) > c, a(x) — \/3(x)\ > c, c < fi(x) < C, then it was stated in [2] that a corresponding 
Feller process exists. Therein was also proposed an example of a mean reverting normal inverse Gaussian- 
like process, a special case of this model with mean is obtained by setting 

a(x) := S(x) := 1, fi(x) := and (3(x) := arctan(x). 
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Note that the mean reversion is not introduced by using simply a drift which drags the process back to 
the origin. It is the choice of j3 which yields an asymmetric distribution that moves the process back to 
the origin. The mean reversion can be observed in Figure |3] which shows samples of the normal inverse 
Gaussian-like process on [0, 1000] with time-step size j^. 

Meixner-like process 

The characteristic function of a Meixner process L t is given by 



where a, r > 0, — ir < b < 71", m E M. Details can be found in (TS] . 

A family of Meixner processes which corresponds to a Feller process can be constructed by substituting 
the parameters a, &, r, m by arbitrary often differcntiablc bounded functions a(x), b(x), r(x), m(x) with 
bounded derivatives. The functions have to be bounded away from the critical values, i.e. < ao < 
a(x), < ro < r(x) and — 7r < b- < b(x) < 6+ < tt for some fixed tq, ao, 6— , &+. For further details 
see [6]. 

Figure H] shows some samples of the Meixner-like process on [0, 100] with time-step size jL and 



The chosen functions satisfy the existence conditions from above. Furthermore the function x t— > a(x) 
yields that the process moves with bigger steps around the origin, to be precise: the Meixner distribution 
has semiheavy tails |17j and the parameter a determines the rate of the exponential decay factor for the 
density. The effect on the sample path can be observed in Figure 0J 

Conclusion 

Using the presented mixture approach one can easily construct Feller models based on given Levy models. 
In these cases the existence of the process is granted. 

Furthermore the presented approximation is a very intuitive way to generate the sample path of a Feller 
processes. Obviously the method requires that one can simulate the increments of the corresponding Levy 
processes. But for Levy processes used in applications, especially together with Monte Carlo techniques, 
this poses no new restriction. 

Thus all necessary tools are available to use Feller processes as models for a wide range of applications. 

Materials and Methods 

The simulations where done in R [55] and the source code of the figures can be found as supporting 
information (Appendix SI). 
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Figure 1. Brownian-Poisson-Cauchy-mixture Feller process. Around the origin it behaves like 
a Poisson process, for smaller values like Brownian motion and a for larger values like a Cauchy process. 




Figure 2. Stable-like processes with a(x) := 1 + y§ ((j - [_3_l) A (Til - f)) - Each position is 
color coded by the corresponding value of the exponent. In the yellow regions it behaves almost like 
Brownian motion, in the red regions it behaves almost like a Cauchy process. 
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Figure 3. Normal inverse Gaussian-like processes with /3(x) 
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Figure 4. Meixner-like process with a(x) := 1 + lOe 25-^ l(_ 5 5 -)(x). The process moves with 
bigger steps around the origin than for larger (and smaller) values. In fact by the choice of x 1— > a(x) the 
rate of the exponential decay of the transition density is reduced around the origin. 
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Appendix SI 

################################################################################ 

# Simulations of Feller processes - B. Boettcher 

################################################################################ 

# Code by B. Boettcher (bjoern. boettcher {at} tu-dresden . de) for the paper: 

# Feller processes: The Next Generation in Modeling. Brownian Motion, Levy 

# processes and beyond. 
# 

# Code executed within R (http://www.r-project.org/), version 2.11.1 
# 

# Note: We opted for clarity, rather than for the most efficient code. The 

# functions are written for an educated user, i.e. the functions do in general 

# not check if the entered parameters have appropriate values . 
# 

# If you want to use this code to simulate Feller processes - instead of just 

# reproducing the figures of the paper - please note: 
# 

# rf ellerprocess (maxtime , timestepsize , startpoint, rldstep) generates the 

# sample path and rldstep is one of 

# rlevymixture - for a process given by the mixture approach 

# rstablef amily - for a stable-like process 

# rnigfamily - for a normal inverse Gaussian-like process 

# rmeixnerf amily - for a Meixner-like process 
# 

# The components of the mixture can be modified at the beginning of the section 

# "the Levy mixture approach" . 

# The families can be modified in the functions rstablef amily , rnigfamily and 

# rmeixnerf amily , respectively. 
# 

# Examples : 

# plot (rf ellerprocess (10, . 001 , , rlevymixture) ) 

# plot(rf ellerprocess (10, 0.001,0, rstablef amily)) 

# plot (rf ellerprocess (1000 ,0 . 1 ,0 ,rnigf amily) ) 

# plot (rf ellerprocess (1000 ,0 . 1 , , rmeixnerf amily) ) 
# 

################################################################################ 

rf ellerprocess <- function (maxtime , timestepsize , startpoint .rldstep) { 

# generates the approximation to the path of the Feller process 
# 

# maxtime = time up to which the path will be simulated 

# timestepsize = time-step size of the increments 

# startpoint = starting point of the process 

# rldstep = a one dimensional function for one step generation 

# (the function should accept two parameters: 

# the time-step size and the current position) 
# 

# returns a list containing the time and position of the process as xy. coords 
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n = ceiling(maxtime/timestepsize) 
v = numeric (n+1) 
v[l] = startpoint 
for ( i in l:n) { 

v[i+l] = rldstep(timestepsize,v[i] ) + v[i] 

} 

invisible (xy . coords ( (0 :n) * timestepsize , v,xlab="time" ,ylab=" space") ) 



################################################################################ 

# the Levy mixture approach 

################################################################################ 

# (here we work with global variables for convenience) 
# 



n = 3 

rlevy = list(n) 
a = list(n) 

rlevy [[1]] = function(t) rnorm(l , sd=sqrt (t) ) 

rlevy [[2]] = function(t) rpois (1 , lambda=t) 

rlevy [[3]] = function(t) rcauchy (1 , scale=t) 



# Number of processes to mix 

# the processes 

# the coefficients 

# L~l = Brownian motion 

# L~2 = Poisson process 

# L~3 = Cauchy process 



# Now one could define just the functions directly (and the process exists 

# always if the functions are non-negative, bounded and Lipschitz continuous) 

# a[[l]] = 

# a[[2]] = 
# 

# Instead, we define regions and generate the functions by linear interpolation 



regions = c( 

-6, # right endpoint of the region where the process is just L~l 

-4,4, # left & right endpoint of the region where the process is just L~2 

6 # left endpoint of the region where the process is just L~3 

) 



# The following yields functions 'a' such that they are =1 on their region and 

# =0 on the others. Constructed by linear interpolation. 

a[[l]] = approxfun(c (regions [1] , regions [2] ), c(l ,0) , method= "linear" , 1 ,0) 
ford in 2: (n-1)) { 

a[ [i] ] =approxfun(c (regions [i*2-3] .regions [i*2-2] , regions [i*2-l] , regions [i*2] ) 
, c (0,1, 1,0) ,method=" linear" ,0,0) 

} 

a[[n]] = approxfun(c (regions [n*2-3] .regions [n*2-2] ) ,c(0, 1) , method= "linear" ,0, 1) 
rlevymixture = function(t,x) { 

# generates an increment of the Levy family constructed as mixture with the 

# above defined processes and coefficients 
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# 

# t = time-step size 

# x = current position 

y = o 

for(i in l:n) 

y = y + a[[i]] (x) *rlevy [ [i] ] (t) 
return(y) 

} 



################################################################################ 

# the symmetric-alpha-stable family 

################################################################################ 
rstable = function(n,t,alpha=l) ■[ 

# generation of increments of an alpha stable process 
# 

# n = number of increments 

# t = time-step size 

# alpha = index of stability 
# 

# We use 

# * the algorithm form Chambers, J. M . ; Mallows, C. L. & Stuck, B. W. 

# A Method for Simulating Stable Random Variables, Journal of the American 

# Statistical Association, 1976, 71, 340-344 

# * the stability, i.e.: 

# X has char.exp. |xi|~alpha ==> t~(l/alpha) X has char.exp. t |xi|~alpha 
W = rexp(n) 

Theta = runif (n, -pi/2, pi/2) 

return (t~ (1/ alpha)* (sin(alpha*Theta) /cos (Theta) " (1/alpha) 
* (cos ( ( 1-alpha) *Theta) /W) ~ ( (1-alpha) /alpha) ) ) 

} 

rstablef amily = function (t,x){ 

# generates an increment of the symmetric alpha stable process family 
# 

# t = time-step size of the increment 

# x = current position 

exponent = function(x) l+19/10*pmin( (x/4-f loor (x/4) ) , (ceiling(x/4) -x/4) ) 

# < exponent < 2 

rstable (1 ,t , exponent (x) ) 

} 



################################################################################ 
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# the normal inverse Gaussian family 

################################################################################ 

rinversegaussian = f unction(t , a,b) { 

# generates an increment of the inverse Gaussian process 
# 

# t = time-step size 

# a,b = parameters of the inverse Gaussian process 
# 

# Algorithm from: V. Seshadri, The inverse Gaussian distribution: a case study 

# in exponential families, Oxford University Press, 1993 (page 203) 



y = rnorm(l)~2 
u = runif(l) 
m = t*a 



xtemp = m+nT2*y/(2*b)- m*sqrt(4*m*b*y+m~2*y~2)/(2*b) 
if (u <= m/(m+xtemp)) return (xtemp) 
else return (m~2/xtemp) 

> 

rnormalinversegaussian = function(t, alpha, beta, delta) { 

# generates an increment of the normal inverse Gaussian process 
# 

# t = time-step size 

# alpha, beta, delta = parameters of the normal inverse Gaussian process 
# 

# Simulation based on subordination of Brownian motion by an inverse Gaussian 

# process 

y = rinversegaussian (t , delta, sqrt (alpha~2-beta~2) ) 
return (bet a*y+rnorm ( 1 , sd=sqrt (y) ) ) 

> 

rnigfamily = function(t ,x) { 

# generates an increment of the normal inverse Gaussian family 
# 

# t = time-step size 

# x = current position 

alpha = function(x) 1 # alpha > 

beta = function(x) -l/pi*atan(x) # -alpha < beta < alpha 

delta = function(x) 1 # delta > 

mu = function(x) # mu is a real number 



return (mu(x) *t+rnormalinversegaussian(t , alpha (x) ,beta(x) ,delta(x) ) ) 

} 
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################################################################################ 

# the Meixner family 

################################################################################ 

# we use the package cxxPack for an implementation of the complex valued 

# Gamma function "cgammaO" 

# to install the package use: 

# install. packages ("cxxPack") 
# 

# Please note, that we have no relation to the author of that package, 
library ( " cxxPack" ) 

rmeixner = function(t,a,b=0,m,r) { 

# generates an increment of the Meixner process 
# 

# t = time-step size 

# a,b,m,r = parameters of the Meixner process 
# 

# The code is based on the method proposed in the preprint : 

# Reiichiro Kawai , Parameter Sensitivity Estimation for Meixner Distribution and 

# Levy Processes, University of Leicester 
# 

# We only implemented the method for b=0 and checked the results of the proposed 

# method by an optical comparison of the resulting histogram with the known 

# density function for various parameters: 
# 

# dmeixner = f unction(t ,a,b,m,r ,x) { 

# density function of the increments of the Meixner process 
# 

# r=r*t 

# m=m*t 

# return ( (2*cos (b/2) ) " (2*r) / (2*a*pi*gamma(2*r) ) 

# *exp (b* (x-m) / a) *abs (cgamma(r+li* (x-m) /a) )~2) 

# } 
# 

# t = 

# a = 

# m = 

# r = 
# 

# x=replicate( 10000, rmeixner (t ,a,0,m,r) ) 

# hist(x,f req=F) 

# f = function(x) dmeixner(t,a,0,m,r,x) 

# curve(f ,col="red" ,add=T) 
# 



if (b!=0) { 

print ( "Error : Only the case b=0 is implemented") 
return ( "Error") 

> 



20 



r=t*r 
m=t*m 
repeat { 

Q = runif (1 ,min=-l ,max=l)/runif (1 ,min=-l ,max=l) 

V = a/2*max(sqrt(2*r) ,2*r)*Q 

U = runif (1) 

if (abs(Q) < 1) { 

if (gamma(r)~2*U < abs (cgamma(r+li*V/a) ) ~2) { 
break 

} 

} 

else { 

if ( ! is . na(cgamma(r+li*V/a) ) ) 

# Note: if abs(V) is to large then cgamma yields NaN (but in fact the value 

# is very small, so also the following inequality would be false. So we 

# just treated the NaN exception separately.) 

if (max(l,2*r)*a~2*gamma(r+l)~2*U/(2*r) < abs (cgamma(r+li*V/a) ) ~2 *V~2){ 
break 

} 

} 

} 

return (V+m) 

> 

rmeixnerf amily = function(t ,x) { 

# generates an increment of the Meixner family 
# 

# t = time-step size 

# x = current position 

a = function(x) if (abs(x)>=5) return(l) else return(l+10*exp(-l/(25-x~2) ) ) 
# a>0 

b = function(x) # -pi < b < pi (implemented only for b=0 ! !) 

m = function(x) # m is a real number 

r = function(x) 1 # b>0 

return(rmeixner(t,a(x) ,b(x) ,m(x) ,r(x) ) ) 

################################################################################ 
################################################################################ 

# Generation of the Figures 

################################################################################ 

# initialization of the random number generator and some graphic parameters 
set.seed(1034) 



par(oma = c (0,0, 0,0), mar=c(4,3.2,2,0.8)) 
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layout (matrix (c(l,2,3,4) , nr ow=2 , byrow=TRUE) , respect=TRUE) 

plottype = "s" 
plotpch = 19 

figplot = function(xy , . . . ) { 

plot(xy,xlab="" ,ylab="" ,pty=plottype,pch=plotpch, . . .) 
title (xlab="time" ,ylab=" space" ,line=2 . 2) 

} 

# the levy mixture 

figplot (rfellerprocess (20 , 0.01,0, rlevymixture) ,main=expression(x [0] == 0)) 
figplot (rfellerprocess (20 , 0.01,0, rlevymixture) ,main=expression(x [0] == 0)) 
figplot (rfellerprocess (20 ,0 . 01 , -10 .rlevymixture) ,main=expression(x [0] == -10) ) 
figplot (rfellerprocess (20 ,0 . 01 , 10, rlevymixture) ,main=expression(x [0] == 10) ) 
dev. copy2eps(f ile="Brownian-Poisson-Cauchy-mixture . eps" , width = 4.8, height = 4.8) 

# the stable-like process 

layout (matrix (c (1 ,2,3,4,5,5) ,nrow=3,byrow=TRUE) , height s=c (1 ,1,0.4) , respect=TRUE) 
par(oma = c (0,0, 0,0), mar=c(4,3.2,2,0.8) ,cex=0.83) 

colorexponent = function(x) pmin( (x/4-f loor (x/4) ) , (ceiling(x/4) -x/4) ) 
cols = heat . colors (100) 

path=rf ellerprocess (20 ,0 . 01 ,0, rstablef amily) 
figplot (path, main=expression(x [0] == 0) 

, col=cols [l+ceiling( colorexponent (path$y) *150)] ) 
path=rf ellerprocess (20 ,0.01,0, rstablef amily) 
figplot (path, main=expression(x [0] == 0) 

, col=cols [l+ceiling( colorexponent (path$y) *150)] ) 
path=rf ellerprocess (20 ,0 . 01 , -2 , rstablef amily) 
figplot (path, main=expression(x [0] == -2) 

, col=cols [l+ceiling( colorexponent (path$y) *150)] ) 
path=rf ellerprocess (20 ,0 . 01 , 2, rstablef amily) 
figplot (path, main=expression(x [0] == 2) 

, col=cols [l+ceiling(colorexponent (path$y) *150)] ) 

# legend for the index of the stable-like process 
plotlegend = functionO { 

par(mar=c(3,0,2,0)) 

plot (c ( 1 , 2) , c (0 , 1) , type="n" , axes=0 ,xlab=" " ,ylab=" " ) 

title(expression(paste("color coding of the values of ", alpha, "(",x, ")"))) 

axis(l) 

n=80 

xr = seq(l,n)/n+l 
xl = xr-l/n 
yb = rep(0,n) 
yt = rep(l,n) 

rect (xl ,yb,xr ,yt , col=cols [1 :n] ,border= NA) 

} 
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plotlegendO 

dev. copy2eps (f ile="stable-like . eps" , width = 4.8, height = 5.76) 

# the normal inverse Gaussian-like process 
par(oma = c (0,0, 0,0), mar=c (4,3 . 2 ,2 ,0 . 8) ) 

layout (matrix (c(l,2,3,4) , nr ow=2 , byrow=TRUE) , respect=TRUE) 

f igplot(rf ellerprocess(1000,0. 1 ,0,rnigf amily) ,main=expression(x [0] == 0)) 

f igplot(rfellerprocess( 1000,0. 1 ,0,rnigf amily) ,main=expression(x [0] == 0)) 

f igplot (rf ellerprocess (1000 ,0 . 1 , -100 ,rnigf amily) ,main=expression(x [0] == -100) ) 

f igplot (rf ellerprocess (1000 ,0 . 1 , 100 ,rnigf amily) ,main=expression(x [0] == 100) ) 

dev. copy2eps (f ile="normal-inverse-gaussian-like . eps" , width = 4.8, height = 4.8) 

# the Meixner-like process 

f igplot (rf ellerprocess (100 ,0 . 1 ,0 ,rmeixnerf amily) ,main=expression(x [0] == 0) ) 
f igplot (rf ellerprocess (100 ,0 . 1 ,0 ,rmeixnerf amily) ,main=expression(x [0] == 0) ) 
f igplot (rf ellerprocess (100 ,0 . 1 , -10 ,rmeixnerf amily) ,main=expression(x [0] == -10) ) 
f igplot (rf ellerprocess (100 ,0 . 1 , 10, rmeixnerf amily) ,main=expression(x [0] == 10) ) 
dev. copy2eps (f ile="meixner-like . eps" , width = 4.8, height = 4.8) 



